% fields H ans E for rho > a & rho < a

function [Hrho, Hphi, Hz, Erho, Ephi, Ez] = func_HE_outside(r, m, CONSTS, coef, p1)

    k0    = CONSTS.k0;
    a     = CONSTS.a;
    eps_a = CONSTS.eps_a;

%     if (isnan(p1))
%         Hrho = 0;
%         Hphi = 0;
%         Hz   = 0;
%         Erho = 0;
%         Ephi = 0;
%         Ez   = 0;
%         return;
%     end
    
    C1 = coef(3);
    C2 = coef(4);

    s  = func_s_from_p(p1, CONSTS);

    S = k0*a*s;
    besselK_m = besselmx(double('K'), m, S*r, 0); %besselk(m, S*r);
    besselK_mp1 = besselmx(double('K'), m+1, S*r, 0); %besselk(m+1, S*r);
    
    Hrho = 1i*(C1*(m./(S*r)).*besselK_m - ...
                        C2*p1*(besselK_mp1 - (m./(S*r)).*besselK_m));
    
    Hphi = -C1*(besselK_mp1 - (m./(S*r)).*besselK_m) + ...
                     C2*p1*(m./(S*r)).*besselK_m;
    
    Hz   = C2*s*besselK_m;    
    
    
    Erho = -C1*(p1/eps_a)*(besselK_mp1 - (m./(S*r)).*besselK_m) + ...
                     C2*(m./(S*r)).* besselK_m;
    
    Ephi = -1i*C1*(p1/eps_a)*(m./(S*r)).*besselK_m + ...
                     1i*C2*(besselK_mp1 - (m./(S*r)).*besselK_m);
    
    Ez   = -(1i/eps_a)*C1*s*besselK_m;
    
end
